Early migration phenology models

As a reminder, the figure below displays a mile-high view of migration timing across different striped bass spawning groups: Hudson River (HR), Delaware River (DE), and Potomac River (PR). The length of the line displays the range of residence of each group in either the natal estuary (colors) or coastal Massachusetts (yellow). The circle shows the median day of arrival in the given array, and the triangle shows the median date of departure. Only fish that were larger than 800 mm at the time of tagging are shown.

Warning: Invalid .internal.selfref detected and fixed by taking a (shallow)
copy of the data.table so that := can add this new column by reference. At an
earlier point, this data.table has been copied by R (or was created manually
using structure() or similar). Avoid names<- and attr<- which in R currently
(and oddly) may copy the whole data.table. Use set* syntax instead to avoid
copying: ?set, ?setnames and ?setattr. If this message doesn't help, please
report your use case to the data.table issue tracker so the root cause can be
fixed or this message improved.

I’m going to make a quick-and-dirty model to compare the timing of entry into coastal Massachusetts by spawning group. The full model includes day of year (“yday”, numeric), spawning group (“group”, factor), the year (“year”, ordered factor), and all interactions.

The response was modeled as binomial, with the cumulative number of fish that have entered Massachusetts waters per group-year combination considered “successes” and the total number of fish that arrived in coastal MA per group-year combination as the number of trials.

Using step-wise variable selection…

Start:  AIC=876.21
cbind(success, trials - success) ~ 0 + group + yday + year + 
    group:yday + group:year + year:yday + group:year:yday

                  Df Deviance    AIC
<none>                 204.73 876.21
- group:yday:year  6   220.74 880.23

…the variables selected are: group, yday, year, group:yday, group:year, yday:year, group:yday:year.

Call:
glm(formula = cbind(success, trials - success) ~ 0 + group + 
    yday + year + group:yday + group:year + year:yday + group:year:yday, 
    family = "binomial", data = binom_data[array == "MA Coast"])

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.9854  -0.6261   0.0378   0.4881   2.1411  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
groupDE             -15.875824   1.328794 -11.948  < 2e-16 ***
groupHR             -11.264119   1.173194  -9.601  < 2e-16 ***
groupPR             -18.541642   2.964276  -6.255 3.97e-10 ***
yday                  0.097072   0.008390  11.570  < 2e-16 ***
year.L               -8.841642   2.298368  -3.847  0.00012 ***
year.Q               -4.173475   2.657588  -1.570  0.11632    
year.C                4.614116   2.973727   1.552  0.12075    
groupHR:yday         -0.033435   0.010911  -3.064  0.00218 ** 
groupPR:yday          0.014470   0.019871   0.728  0.46650    
groupHR:year.L        8.480084   3.853682   2.201  0.02777 *  
groupPR:year.L        9.320778   7.243423   1.287  0.19817    
groupHR:year.Q       -2.710364   3.545181  -0.765  0.44456    
groupPR:year.Q        0.729079   6.496961   0.112  0.91065    
groupHR:year.C       -4.344584   3.207141  -1.355  0.17553    
groupPR:year.C        4.656374   5.652772   0.824  0.41009    
yday:year.L           0.061693   0.014208   4.342 1.41e-05 ***
yday:year.Q           0.023110   0.016779   1.377  0.16842    
yday:year.C          -0.033881   0.019006  -1.783  0.07464 .  
groupHR:yday:year.L  -0.055324   0.023270  -2.377  0.01743 *  
groupPR:yday:year.L  -0.072239   0.044520  -1.623  0.10467    
groupHR:yday:year.Q   0.018236   0.021822   0.836  0.40333    
groupPR:yday:year.Q   0.002267   0.039743   0.057  0.95451    
groupHR:yday:year.C   0.033170   0.020270   1.636  0.10177    
groupPR:yday:year.C  -0.013135   0.034307  -0.383  0.70181    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2656.34  on 202  degrees of freedom
Residual deviance:  204.73  on 178  degrees of freedom
AIC: 876.21

Number of Fisher Scoring iterations: 5

When looking at pairwise contrasts by year, the only significant differences occur in 2018, where arrival in MA was more rapid in DE and PR fish when compared to HR fish (greater yday slope).

year = 2016:
 contrast estimate      SE  df z.ratio p.value
 DE - HR  -0.00538 0.00971 Inf  -0.554  0.8443
 DE - PR  -0.06700 0.05253 Inf  -1.275  0.4091
 HR - PR  -0.06162 0.05187 Inf  -1.188  0.4604

year = 2017:
 contrast estimate      SE  df z.ratio p.value
 DE - HR   0.00793 0.00513 Inf   1.547  0.2691
 DE - PR  -0.02068 0.02399 Inf  -0.862  0.6644
 HR - PR  -0.02861 0.02372 Inf  -1.206  0.4494

year = 2018:
 contrast estimate      SE  df z.ratio p.value
 DE - HR   0.07717 0.02762 Inf   2.794  0.0144
 DE - PR  -0.00599 0.03971 Inf  -0.151  0.9875
 HR - PR  -0.08317 0.02924 Inf  -2.845  0.0124

year = 2019:
 contrast estimate      SE  df z.ratio p.value
 DE - HR   0.05401 0.03196 Inf   1.690  0.2090
 DE - PR   0.03579 0.03750 Inf   0.955  0.6056
 HR - PR  -0.01822 0.04316 Inf  -0.422  0.9064

P value adjustment: tukey method for comparing a family of 3 estimates 

NOTE: Results may be misleading due to involvement in interactions

$emmeans
year = 2016:
 group emmean     SE  df asymp.LCL asymp.UCL
 DE     1.376 0.2072 Inf     0.970     1.782
 HR     0.953 0.0709 Inf     0.814     1.092
 PR     4.711 1.7653 Inf     1.251     8.171

year = 2017:
 group emmean     SE  df asymp.LCL asymp.UCL
 DE     0.738 0.1255 Inf     0.492     0.984
 HR     0.364 0.0666 Inf     0.233     0.495
 PR     2.738 1.0116 Inf     0.755     4.721

year = 2018:
 group emmean     SE  df asymp.LCL asymp.UCL
 DE     4.691 1.0661 Inf     2.602     6.780
 HR     0.575 0.1176 Inf     0.344     0.805
 PR     1.778 0.4627 Inf     0.871     2.685

year = 2019:
 group emmean     SE  df asymp.LCL asymp.UCL
 DE     4.627 0.5827 Inf     3.485     5.769
 HR     2.177 0.7645 Inf     0.678     3.675
 PR     2.711 0.9279 Inf     0.893     4.530

Results are averaged over the levels of: yday 
Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 

$contrasts
year = 2016:
 contrast estimate    SE  df z.ratio p.value
 DE - HR     0.423 0.219 Inf   1.929  0.1304
 DE - PR    -3.335 1.777 Inf  -1.876  0.1455
 HR - PR    -3.758 1.767 Inf  -2.127  0.0844

year = 2017:
 contrast estimate    SE  df z.ratio p.value
 DE - HR     0.374 0.142 Inf   2.635  0.0229
 DE - PR    -2.000 1.019 Inf  -1.962  0.1218
 HR - PR    -2.374 1.014 Inf  -2.342  0.0502

year = 2018:
 contrast estimate    SE  df z.ratio p.value
 DE - HR     4.116 1.073 Inf   3.838  0.0004
 DE - PR     2.913 1.162 Inf   2.506  0.0327
 HR - PR    -1.203 0.477 Inf  -2.521  0.0314

year = 2019:
 contrast estimate    SE  df z.ratio p.value
 DE - HR     2.450 0.961 Inf   2.549  0.0291
 DE - PR     1.915 1.096 Inf   1.748  0.1875
 HR - PR    -0.535 1.202 Inf  -0.445  0.8967

Results are averaged over the levels of: yday 
Results are given on the log odds ratio (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

Entry into MA waters significantly different between DE and HR fish, but not between PR and HR or PR and DE fish.

Start:  AIC=991.5
cbind(success, trials - success) ~ 0 + group + yday + group:yday + 
    year

             Df Deviance     AIC
<none>            382.39  991.50
- group:yday  2   394.16  999.26
- year        3   528.96 1132.06

Call:
glm(formula = cbind(success, trials - success) ~ 0 + group + 
    yday + group:yday + year, family = "binomial", data = binom_data[array == 
    "MA Coast"])

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0910  -0.9008   0.0119   0.9630   5.2940  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
groupDE       -8.542018   0.499708 -17.094  < 2e-16 ***
groupHR      -11.015699   0.352014 -31.293  < 2e-16 ***
groupPR      -10.413725   1.745313  -5.967 2.42e-09 ***
yday           0.033748   0.001915  17.623  < 2e-16 ***
year2017      -0.374309   0.070063  -5.342 9.17e-08 ***
year2018       0.906559   0.108262   8.374  < 2e-16 ***
year2019       0.300201   0.130164   2.306 0.021092 *  
groupHR:yday   0.008022   0.002301   3.486 0.000491 ***
groupPR:yday   0.007858   0.007130   1.102 0.270393    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2536.41  on 192  degrees of freedom
Residual deviance:  382.39  on 183  degrees of freedom
AIC: 991.5

Number of Fisher Scoring iterations: 5

$emtrends
 group yday.trend      SE  df asymp.LCL asymp.UCL
 DE        0.0337 0.00191 Inf    0.0300    0.0375
 HR        0.0418 0.00133 Inf    0.0392    0.0444
 PR        0.0416 0.00689 Inf    0.0281    0.0551

Results are averaged over the levels of: year 
Confidence level used: 0.95 

$contrasts
 contrast  estimate      SE  df z.ratio p.value
 DE - HR  -0.008022 0.00230 Inf  -3.486  0.0014
 DE - PR  -0.007858 0.00713 Inf  -1.102  0.5127
 HR - PR   0.000164 0.00700 Inf   0.023  0.9997

Results are averaged over the levels of: year 
P value adjustment: tukey method for comparing a family of 3 estimates 

exit times not significantly different